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Abstract 

ICM-Lite  is  envisioned  as  a  tool  to  rapidly  predict  and  quantify  ecosystem 
benefits  in  aquatic  systems.  ICM-Lite  consists  of  three  modules:  a  water 
quality  module,  a  graphical  user  interface,  and  an  ecosystems  benefits 
module.  This  publication  documents  the  formulations  of  the  water  quality 
module.  Documentation  includes  the  mass  balance  equation,  kinetics  in 
the  water  column,  and  representations  of  sediment-water  fluxes.  The  mass 
balance  relationships  are  intended  to  incorporate  flows  and  volumes 
provided  by  the  user  and  obtained  outside  of  ICM-Lite.  The  kinetics 
include  representations  of  the  aquatic  carbon,  nitrogen,  phosphorus,  and 
oxygen  cycles.  Salinity,  temperature,  suspended  solids,  and  one  user- 
defined  substances  are  included,  as  well.  Sediment-water  fluxes  of  organic 
matter,  ammonium,  nitrate,  phosphate,  and  dissolved  oxygen  are 
considered. 
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1  Introduction 

1.1  The  problem 

The  U.S.  Army  Corps  of  Engineers  (USAGE)  needs  to  quantify  ecosystem 
benefits  to  organisms  and  to  the  aquatic  environment  that  result  from 
USAGE  restoration  projects.  The  USAGE  presently  operates  a  number  of 
predictive  hydrodynamic  and  water  quality  models  which  compute 
environmental  properties  such  as  water  temperature,  water  clarity,  and 
dissolved  oxygen.  While  these  tools  are  superb,  their  applicability  can  be 
limited  in  the  present  operational  climate.  One  limitation  is  their 
complexity.  The  existing  models  provide  results  at  high  levels  of  spatial 
and  temporal  resolution  by  employing  detailed  representations  of  physical 
and  biogeochemical  processes.  The  production  of  results  at  these  levels  of 
detail  requires  large  volumes  of  data,  for  setup  and  validation,  and 
extensive  time  for  application  by  skilled  specialists.  These  resource 
demands  result  in  corresponding  high  costs.  Both  the  time  and  the  cost 
weigh  against  the  current  trends  in  planning  as  expressed,  for  example,  in 
the  recent  “3x3x3”  directive  (Walsh  2012). 

Resource  requirements,  especially  time,  comprise  a  second  limitation  to 
the  utility  of  current  planning  models.  They  are  ill  suited  for  uncertainty 
analysis.  Multiple  approaches  to  quantifying  uncertainty  associated  with 
ecosystem  benefits  are  available  but  all  share  in  common  the  need  to 
execute  and  analyze  large  numbers  of  model  runs  based  on  varying  input 
parameter  sets.  The  time  involved  to  set  up,  execute,  and  analyze  results 
from  the  current  suite  of  predictive  tools  precludes  their  employment  in 
Monte  Garlo  analyses  and  other  examinations  of  uncertainty. 

An  additional  limitation  in  existing  planning  tools  is  inherent  in  the  nature 
of  the  information  produced.  The  current  models  provide  basic  physical 
properties,  such  as  temperature  or  turbidity,  at  locations  and  times 
determined  by  model  time-step  and  computational  grid.  The  information 
desired,  however,  is  often  characterized  in  terms  of  ecosystem  benefits  or 
habitat  suitability  expressed  in  units  and  over  spatial  and  temporal  scales 
that  are  not  explicit  in  the  model  computational  details. 
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The  USAGE  will  always  have  projects  that  demand  the  state-of-the-art  in 
high-fidelity  planning  models.  Clearly,  however,  the  USAGE  also  requires  a 
simpler  tool  (or  tools)  with  the  following  characteristics: 

•  fast  and  easy  to  apply 

•  capable  of  uncertainty  analysis 

•  applicable  at  the  District  level  by  nonspecialists 

•  results  provided  in  useful  and  immediately  understandable  terms. 

1.2  This  study 

The  present  study  consists  of  three  phases: 

•  reduce  the  complexity  of  a  popular  USAGE  planning  tool  to  rapidly 
provide  basic  results 

•  wrap  the  simplified  tool  in  a  graphical  interface  (GUI)  to  facilitate  use 
at  the  District  level 

•  enable  the  GUI  to  conduct  uncertainty  analyses  and  to  provide  output 
in  desired  terms. 

The  USAGE  ERDC  is  the  creator  and  operator  of  the  Corps  of  Engineers 
Water  Quality  Integrated  Compartment  Model  (CE-QUAL-ICM  or  simply 
ICM).  This  state-of-the-art  model  has  been  employed  in  lakes,  rivers,  and 
estuaries  throughout  the  United  States  and  overseas.  The  model  is  subject 
to  the  limitations  noted  above,  however.  ICM  requires  experienced  users 
and  can  be  expensive  and  time  consuming  to  apply.  The  initial  phase  of 
this  project  will  create  an  ICM-Lite,  a  simplified  version  of  ICM  kinetics. 
ICM-Lite  will  compute  aquatic  properties  including  salinity,  temperature, 
dissolved  oxygen,  chlorophyll,  total  suspended  solids,  and  turbidity.  ICM- 
Lite  is  envisioned  for  an  aquatic  system  that  is  represented  by  one  well- 
mixed  cell  or  by  a  simple  network  of  cells.  Basic  hydrodynamics  will  be 
based  on  flows  and  exchange  coefficients. 

The  second  phase  of  the  project  will  be  the  development  of  a  GUI  that  will 
facilitate  ICM-Lite  application  by  the  novice  user.  Both  GUI  and  model 
will  operate  on  a  Windows  PC.  The  GUI  is  envisioned  to  prompt  the  user 
for  desired  information  and  to  suggest  default  values  when  these  are 
appropriate.  Habitat  criteria  or  similar  information  input  to  the  GUI  by 
the  user  will  be  employed  to  translate  the  native  information  produced  by 
the  model,  such  as  dissolved  oxygen,  into  ecosystem  benefits  or  other 
useful  terms. 
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The  swift  execution  of  ICM-LITE  will  allow  for  Monte  Carlo  simulations 
and  uncertainty  analysis.  The  GUI  will  provide  built-in  guidance  for  this 
process.  As  a  result,  ecosystem  benefits  and/or  restoration  of  desirable 
habitat  can  be  expressed  as  an  expected  range  of  values  along  with 
associated  probabilities. 

1.3  This  report 

The  project  plan  calls  for  the  first  project  year  to  be  devoted  to  reducing 
the  existing  CE-QUAL-ICM  from  its  present  30+  state  variables  to  a  suite 
of  8  to  10,  including  basic  physical  properties  such  as  temperature  and 
salinity,  nutrients,  dissolved  oxygen,  and  chlorophyll.  The  present  report 
details  the  formulation  of  the  new  ICM-Lite. 
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2  The  Conservation  of  Mass  Equation 

2.1  The  model  grid 

Application  of  the  model  requires  division  of  the  study  system  into  a  grid 
of  discrete  volumes  or  cells.  Although  each  volume  is  three-dimensional, 
the  grid  may  be  one-,  two-,  or  three-dimensional  (iD,  2D,  3D),  depending 
on  the  arrangement  of  the  cells.  An  example  of  a  2D  grid  is  shown  as 
Figure  1.  This  grid  contains  ten  cells  in  the  longitudinal  dimension,  one 
cell  in  the  lateral  dimension,  and  three  cells  in  the  vertical  dimension. 


Figure  1.  Example  of  2D  grid  (elevation).  The  grid  is  ten  cells  in  the  longitudinal  direction, 
three  cells  in  the  vertical  direction,  and  one  cell  in  the  lateral  direction. 


Inflow 

2 

3 

4 

5 

6 

7 

8 

9 

Tidal 

Boundary 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20  ^ 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30  < 

Each  cell  in  the  grid  is  assigned  a  unique  number  or  index  (Figure  2). 
Interfaces  are  numbered  where  flows  pass  between  cells  or  where  cells 
adjoin  open  boundaries.  Faces  adjacent  to  solid  boundaries  are  not 
numbered.  The  grid  is  unstructured.  That  is,  the  cell  index  contains  no 
information  that  indicates  cell  location  in  a  3D  coordinate  system.  Neither 
is  there  a  general  relationship  between  the  indices  of  adjacent  cells  or 
between  cells  and  flow  faces.  A  connectivity  or  map  file  is  required  that 
locates  cells  and  faces  relative  to  each  other. 

The  unstructured  grid  of  discrete  volumes  provides  maximum  flexibility  in 
grid  configuration.  No  restriction  is  placed  on  cell  shape  or  number  of  flow 
faces  per  cell.  In  the  complete  ICM,  the  flows  between  cells  usually  come 
from  an  accompanying  hydrodynamic  model.  The  map  file  is  constructed 
by  software  appended  to  the  hydrodynamic  model.  For  ICM-Lite,  the  user 
is  anticipated  to  create  the  mapping  for  a  relatively  simple  grid 
configuration.  Flows  are,  likewise,  anticipated  to  be  entered  by  the  user 
into  a  hydrodynamic  file.  These  flows  may  be  extracted  from  a 
hydrodynamic  model  or  derived  otherwise,  for  example,  from  stream 
gauges.  The  configuration  of  the  map  and  hydrodynamic  files  will  be 
described  in  the  users’  guide  to  be  produced  in  the  second  study  year. 
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Figure  2.  The  numbering  system  for  ceiis,  faces,  and  flows.  The  numbering  scheme 
corresponds  to  the  upper  left  cells  in  the  grid  shown  in  Figure  1. 


2.2  The  conservation  of  mass  equation 

The  foundation  of  ICM  is  the  3D  mass  conservation  equation  for  a  control 
volume.  ICM  solves,  for  each  volume  and  for  each  state  variable,  the 
following  equation: 


8^ 

8f 


8C 

8x. 


(1) 


in  which: 

Vj  =  volume  of  jth  control  volume  (m3) 

Cj  =  concentration  in  jth  control  volume  (g  m-3) 

Qk  =  volumetric  flow  across  flow  face  k  of  jth  control  volume 
(m3  s“i) 

Ck  =  concentration  in  flow  across  flow  face  k  (g  m-3) 

Ak  =  area  of  flow  face  k  (m^) 
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Dk  =  diffusion  coefficient  at  flow  face  k  (m^  s  ^) 
n  =  number  of  flow  faces  attached  to  jth  control  volume 
Sj  =  external  loads  and  kinetic  sources  and  sinks  in  ith  control 
volume  (g  S'l) 

t,x=  temporal  and  spatial  coordinates 

2.3  Discretization  of  the  mass  conservation  equation 

Solution  of  the  conservation-of-mass  equation  on  a  digital  computer 
requires  specification  of  parameter  values  and  discretization  of  the 
continuous  derivatives.  Numerous  formulae  for  evaluation  and 
discretization  exist.  Formulae  employed  in  ICM  were  selected  based  on 
computational  efficiency  and  accuracy. 


The  conservation-of-mass  equation  is  solved  in  two  steps.  In  the  first  step, 
an  intermediate  value  is  computed.  The  intermediate  value  includes  the 
effects  of  change  in  cell  volume,  longitudinal  and  lateral  transport,  and 
external  loading.  In  the  second  step,  the  effects  of  vertical  transport  are 
computed. 


2.3.1  Longitudinal  and  lateral  advection 

Solution  to  the  conservation  of  mass  equation  in  the  longitudinal  and 
lateral  directions  is  via  explicit  time-stepping.  That  is: 


C 


nhf  nhf 

E<?»c.+EAOi 

fc=l  it=l 


(2) 


in  which: 

Cf  =  concentration  in  jth  control  volume  after  volume  change, 
loading,  longitudinal/lateral  transport  processes 
=  volume  ofjth  control  volume  at  time  f=Af 
Af  =  discrete  time-step 

nhf  =  number  of  longitudinal  and  lateral  flow  faces  attached  to  jth 
control  volume 


Except  where  noted  by  f+Af,  all  terms  in  Equation  2  are  evaluated  at  time  t. 

Solution  of  Equation  2  requires  evaluation  of  the  concentrations  flowing 
across  cell  faces,  Ck.  Two  options  are  provided  within  ICM.  The  first  is 
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backwards  or  upwind  differencing.  In  upwind  differencing,  concentration 
in  the  flow  across  any  face  is  taken  as  concentration  in  the  cell  upstream  of 
the  face  (Figure  3).  Upstream  is  defined  relative  to  direction  of  the  flow. 
Upstream  has  no  relation  to  the  cell  coordinate  system. 


Figure  3.  Derivation  of  CkaX.  ceii  faces  for  QUICKEST  and  UPWIND  schemes. 


A  second  approximation  to  Ck  fits  a  parabola  to  concentration  in  three 
adjacent  cells  (Figure  3).  For  uniform  grid  spacing, 

C.=|[C,.+C,]-i[C„+C,-2C,]  (3) 

The  parabolic  fit  is  part  of  a  numerical  scheme  known  as  QUICKEST 
(Leonard  1979).  While  upwind  differencing  provides  computational 
simplicity,  the  upwind  method  is  less  accurate  and  less  stable  than 
QUICKEST.  The  primary  disadvantage  of  QUICKEST  is  that  the  method 
sometimes  generates  negative  concentrations  when  advecting  sharp 
concentration  gradients.  A  second  disadvantage  is  that  QUICKEST  cannot 
be  implemented  on  highly  irregular  grids  in  which  two  adjacent  upstream 
cells  cannot  be  readily  identified. 
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Detailed  knowledge  of  the  advective  solution  schemes  is  not  necessary  to 
execute  the  model.  The  upwind  and  QUICKEST  approximations  were 
reviewed  to  indicate  the  information  required  by  the  model  to  compute 
advective  transport  in  the  longitudinal  and  lateral  directions.  To  compute 
advective  transport  in  any  cell,  the  model  requires 

•  cell  volume 

•  indices  of  longitudinal  and  lateral  flow  faces  adjoining  the  cell 

•  indices  of  adjoining  and  next-most  adjoining  cells 

•  volumetric  flow  across  the  indexed  flow  faces 

•  length  of  indexed  cells. 

The  required  information  is  provided  in  the  map,  geometry,  and 
hydrodynamics  files.  Formats  of  these  files  are  provided  in  the  users’  guide. 


2.3.2  Longitudinal  and  lateral  dispersion 

Computation  of  longitudinal  and  lateral  dispersion  requires  discrete 
approximation  of  the  continuous  derivative  in  the  dispersion  term  of 
Equation  2.  The  basic  approximation  is 


8C  Cj-C, 
AX 


(4) 


in  which: 

Ax  =  distance  between  centers  of  two  cells 

A  higher-order  correction  to  the  basic  expression  is  computed  when  the 
QUICKEST  scheme  is  employed. 

2.3.3  Vertical  transport 

Solution  to  the  conservation-of-mass  equation  in  the  longitudinal  and 
lateral  directions  is  by  an  explicit  method.  That  is,  all  parameters  in  the 
discretized  equation  are  evaluated  at  time  t  except  the  unknown  C*.  The 
explicit  method  is  suited  for  transport  dominated  by  advection  rather  than 
diffusion  or  dispersion.  In  the  vertical  direction,  diffusion  is  a  significant 
or  dominant  component  of  transport.  Solution  of  vertical  transport  by  an 
explicit  method  requires  a  small  time-step  and  consumes  large  amounts  of 
computer  time.  In  ICM,  solution  to  vertical  transport  is  by  a  partly  or  fully 
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implicit  scheme  that  practically  frees  the  computation  from  stability 
conditions  imposed  by  vertical  transport. 

The  mass-conservation  equation  in  the  vertical  direction  (Figure  4)  can  be 
expressed  as 


(i-e)-E 


V, 


At 


■^t+At 


A-Dk 

vr 


8C^ 

8z 


(5) 


in  which: 

0  = 
nvf  = 

z  = 

implicit  weighting  factor  (0  <  0  <  1) 
number  of  vertical  faces 
vertical  coordinate. 

The  implicit  weighting  factor,  G,  determines  whether  vertical  advection  is 
computed  explicitly  (0  =  o),  implicitly  (0  =  1),  or  is  weighted  between  the 
two  extremes  (o  <  0  <  1).  Computational  stability  is  enhanced  as  0  — >  1,  at 
the  expense  of  increased  numerical  diffusion.  The  value  0  =  0.75  is 
recommended. 


Since  vertical  velocities  are  usually  much  less  than  longitudinal  velocities, 
the  enhanced  accuracy  of  the  QUICKEST  scheme  is  not  necessary. 
Concentrations  at  cell  interfaces,  Ck  and  are  computed  by  linear 
interpolation  between  concentrations  at  the  centers  of  adjoining  cells: 


A 


C^■Az^+CJ■Az■ 
Az.  +  Azj 


(6) 


The  spatial  gradient  in  the  diffusion  term  is  evaluated  by  central  difference 
evaluated  at  time-step  f+Af. 

The  solution  scheme  for  vertical  transport  is  an  implicit  scheme,  which 
means  that  the  equation  for  concentration  in  any  cell  at  time  f+Af  (e.g.. 
Equation  5)  contains  multiple  unknowns.  Computation  of  concentration  in 
any  one  cell  requires  solution  of  a  set  of  simultaneous  equations  for 
concentrations  in  a  column  of  cells  extending  from  water  surface  to 
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Figure  4.  Definition  of  iNFLOW,  OUTFLOW,  and  TiDAL  boundary  conditions.  C&is  a  user- 

specified  boundary  concentration. 


INFLOW 


OUTFLOW 


TIDAL 


bottom.  Details  of  the  solution  scheme  are  not  necessary  to  operate  the 
model.  The  user  must  provide,  however,  the  following  information 
required  to  compute  vertical  transport: 

•  indices  of  all  cells  in  a  column 

•  indices  of  vertical  flow  faces  adjoining  all  cells  in  a  column 

•  volumes  of  all  cells  in  a  column 

•  volumetric  flow  across  the  indexed  flow  faces 

•  diffusion  coefficients  at  indexed  flow  faces 

•  length  of  indexed  cells. 
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The  required  information  is  provided  in  the  map,  geometry,  and 
hydrodynamics  files.  Formats  of  these  files  are  provided  in  the  users’  guide 
(in  preparation). 

2.3.4  Summary  of  numerical  solution  scheme 

The  model  solves  the  conservation-of-mass  equation  through  a  step-by- 
step  procedure: 

1.  Evaluate  internal  sources  and  sinks.  These  include  kinetics 
transformations  and  sediment-water  fluxes.  This  step  provides  a  partial 
computation  of  ZSj  in  Equation  i. 

2.  Add  effects  of  external  loads.  This  step  completes  computation  of  ZSj  in 
Equation  i. 

3.  Specify  longitudinal  and  lateral  advection  and  diffusion  at  all  interfaces. 
This  step  provides  quantities  required  to  solve  Equation  2. 

4.  Compute  concentration  at  time  f+Af  in  all  cells  resulting  from  volume 
changes,  kinetics,  external  loads,  and  longitudinal/lateral  transport.  This 
step  is  the  solution  to  Equation  2.  Eor  iD  or  2D  (longitudinal/lateral) 
systems,  solution  of  the  conservation-of-mass  equation  is  complete  at  this 
point.  Eor  2D  (longitudinal/vertical)  or  3D  systems,  the  result  is  an 
intermediate  solution  prior  to  computation  of  vertical  transport. 

5.  Compute  vertical  transport  from  surface  to  bottom.  Computation  is  by 
columns.  Each  cell  at  the  water  surface  represents  the  top  of  one  column. 

2.4  ICM  time-step 

Temporal  integration  of  the  conservation-of-mass  equation  (1)  is 
accomplished  in  discrete  time-steps  Af  (Equations  2,  5).  Integration  in 
discrete  steps  provides  an  approximation  to  the  continuous  solution  of  the 
original  differential  equation.  As  Af  — >  o,  the  solution  of  the  approximate 
equation  converges  on  the  solution  of  the  continuous  equation,  although  at 
great  cost  in  computation  time.  As  Af  — >  00,  computation  time  diminishes, 
but  the  solution  of  the  discrete  equation  diverges  from  solution  of  the 
continuous  equation.  Eor  sufficiently  large  Af,  the  numerical  solution  may 
exhibit  large  oscillations  or  instabilities  that  produce  computational 
failures.  The  occurrence  of  instabilities  is  prevalent  in  explicit  rather  than 
implicit  solution  schemes.  Typical  practice  in  numerical  modeling  is  to 
select  the  largest  time-step  possible,  to  minimize  computation  time,  while 
remaining  in  predefined  stability  limits. 
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The  time-step  employed  is  determined  by  an  autostepping  algorithm.  The 
algorithm  computes  permissible  time-step  based  on  flow,  dispersion,  and 
cell  dimension.  As  a  consequence  of  autostepping,  the  time-step  varies 
throughout  a  model  run.  The  time-step  is  always  near  the  maximum 
permissible  time-step.  Autostepping  minimizes  computation  time  while 
meeting  stability  requirements.  For  typical  ICM  applications,  the  time- 
step  is  on  the  order  of  minutes. 

2.5  Boundary  conditions 

Boundary  conditions  must  be  specified  at  the  flow  faces  along  the  edges  of 
the  grid.  Through  these  faces,  material  is  exchanged  with  the  environment 
outside  the  model  domain.  Boundary  flow  faces  are  allowed  only  at  the 
longitudinal  and  lateral  limits  of  the  grid.  No  flow  is  allowed  through  the 
surface  and  bottom.  Cell  faces  at  the  water  surface  and  bottom  are  not 
indexed,  nor  are  cell  faces  indexed  along  longitudinal  and  lateral  edges  of 
the  grid  through  which  flow  does  not  occur. 

Treatment  of  open  boundary  conditions  requires  selection  of  the 
numerical  scheme  and  specification  of  concentration  in  the  environment 
beyond  the  grid.  Three  types  of  boundary  conditions  are  defined:  inflow 
boundaries,  outflow  boundaries,  and  tidal  boundaries  (Figure  4). 

2.5.1  Inflow  boundary  conditions 

Flow  at  inflow  boundaries  is  unidirectional  into  the  model  domain.  Typical 
inflow  boundaries  include  flow  at  an  estuarine  fall  line  and  river  inflow  at 
the  head  of  a  reservoir.  The  model  employs  upwind  differencing  at  inflow 
boundaries.  Upwind  differencing  occurs  whether  or  not  the  QUICKEST 
scheme  is  specified  for  advection  within  the  interior  of  the  grid.  Upwind 
differencing  ensures  that  the  concentration  of  flow  entering  the  domain  is 
the  specified  boundary  concentration.  If  QUICKEST  were  employed  at  an 
inflow  boundary,  the  three-point  weighting  scheme  would  compute  an 
influence  of  concentration  within  the  system  on  concentration  entering  the 
system.  Dispersion  is  set  to  zero  at  inflow  boundaries  to  ensure  that  the 
mass  entering  the  system  is  determined  exactly  as  the  product  of  flow  and 
concentration. 
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2.5.2  Outflow  boundary  conditions 

Flow  at  outflow  boundaries  is  unidirectional  out  of  the  model  domain. 
Typical  outflow  boundaries  include  flow  over  a  spillway  and  flow  at  the 
downstream  end  of  a  modeled  river  reach.  The  model  employs  upwind 
differencing  at  outflow  boundaries.  Upwind  differencing  occurs  whether  or 
not  the  QUICKEST  scheme  is  specified  for  advection  within  the  interior  of 
the  grid.  Upwind  differencing  ensures  that  the  concentration  of  flow 
leaving  the  domain  is  determined  solely  by  conditions  within  the  domain. 
If  QUICKEST  were  employed  at  an  outflow  boundary,  the  three-point 
weighting  scheme  would  compute  an  influence  of  concentration  outside 
the  system  on  concentration  leaving  the  system.  This  situation  is  clearly 
impossible  at  a  spillway,  for  example.  Dispersion  is  set  to  zero  at  outflow 
boundaries  based  on  the  same  reasoning.  Material  leaving  the  system 
cannot  mix  back  into  the  system. 

2.5.3  Tidal  boundaries 

The  typical  tidal  boundary  is  situated  at  the  mouth  of  an  estuary  or  lagoon. 
Elow  at  tidal  boundaries  is  bidirectional,  depending  on  the  phase  of  the 
tide.  Elow  is  into  the  system  on  the  flood  tide  and  out  of  the  system  on  the 
ebb.  Material  is  free  to  mix  in  both  directions  across  the  interface.  The 
model  employs  upwind  or  QUICKEST  at  tidal  boundaries,  depending  on 
the  scheme  specified  for  the  interior  of  the  system.  Dispersion  is  calculated 
at  tidal  boundaries. 
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3  Kinetics  in  the  Water  Coiumn 

The  Lite  version  of  ICM  incorporates  ii  state  variables  in  the  water 
column: 

•  temperature 

•  salinity 

•  fixed  suspended  solids 

•  phytoplankton  (as  carbon) 

•  organic  carbon 

•  ammonium 

•  nitrate 

•  organic  nitrogen 

•  phosphate 

•  organic  phosphorus 

•  dissolved  oxygen. 

An  additional  user-defined  variable  can  be  implemented  to  provide  basic 
representation  of  a  substance,  such  as  bacteria  or  contaminant,  not  in  the 
standard  model  suite.  Additional  quantities  can  be  derived  from  the  state 
variable  suite.  Diffuse  light  attenuation  is  one  important  derived  variable 
which  is  calculated  based  on  computed  organic  and  inorganic  particulate 
matter. 

The  fate  and  transport  of  all  model  state  variables  are  computed  based  on 
the  mass  conservation  equation  (Equation  i)  presented  in  the  preceding 
chapter.  For  notational  simplicity,  the  transport  terms  are  dropped  from 
the  equation  when  describing  the  kinetics  formulations  which  contribute 
to  the  A  Sj. 

3.1  Temperature 

Computation  of  temperature  employs  a  conservation  of  internal  energy 
equation  that  is  analogous  to  the  conservation  of  mass  equation.  For 
practical  purposes,  the  internal  energy  equation  can  be  written  as  a 
conservation  of  temperature  equation.  The  only  source  or  sink  of 
temperature  considered  is  exchange  with  the  atmosphere.  Atmospheric 
exchange  is  considered  proportional  to  the  temperature  difference 
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between  the  water  surface  and  a  theoretical  equilibrium  temperature 
(Edinger  et  al.  1974): 


KT 


pCpH 


■{Te-T) 


(7) 


in  which: 


T  =  water  temperature  (oC) 

Te  =  equilibrium  temperature  (^C) 

KT  =  heat  exchange  coefficient  (watt  m-2  oC  O 
Cp  =  specific  heat  of  water  (4200  watt  s  kg  ^  oC  ^) 
p  =  density  of  water  (1000  kg  m-3) 

H  =  depth  of  water  column  or  surface  cell  (m). 

3.2  Salinity 


Salinity  is  modeled  by  the  conservation  of  mass  equation  with  no  internal 
sources  or  sinks 


3.3  Fixed  solids 

Fixed  solids  comprise  the  inorganic  particles  suspended  in  the  water 
column.  (Total  Suspended  Solids  is  a  derived  quantity  based  on  computed 
fixed  solids  and  organic  matter).  The  only  internal  source  or  sink  of  fixed 
solids  is  settling: 


o  o 

-^FS  =  -Ws~FS  (8) 

5f  5z 


in  which: 

FS  =  fixed  solids  concentration  (g  m-3) 

Ws  =  fixed  solids  settling  rate  (m  d'O 
z  =  vertical  coordinate. 

3.4  Light  attenuation 

Irradiance  decreases  exponentially  as  a  function  of  depth  below  the  water 
surface: 
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l{z)  =  Io-e-^^'^  (9) 

in  which: 

/(z)  =  irradiance  at  depth  z  (E  m-2  d  O 
lo  =  irradiance  at  water  surface  (E  m-2  d  ^ 

Ke  =  diffuse  light  attenuation  coefficient  (m-i)- 

The  diffuse  light  attenuation  coefficient  is  computed  by  a  partial 
attenuation  model  based  on  computed  suspended  solids: 

Ke  =  al  +  a2FS  +  a3VSS  (10) 

in  which: 

Ke  =  coefficient  of  diffuse  light  attenuation  (m-i) 
ai  =  background  attenuation  (m-i) 
a2  =  attenuation  by  inorganic  suspended  solids  (m^  g  i) 
as  =  attenuation  by  organic  suspended  solids  (m^  g-i) 

FS  =  inorganic  (fixed)  suspended  solids  concentration  (g  m-3) 

VSS  =  organic  (volatile)  suspended  solids  concentration  (g  m-3). 

The  empirical  coefficients  ai-as  are  provided  by  the  user.  VSS  is  based  on 
computed  organic  matter  concentration.  The  units  of  irradiance  and  of  Ke 
are  determined  by  requirements  of  the  phytoplankton  model.  Often,  water 
quality  or  habitat  criteria  are  based  on  turbidity,  which  is  a  measure  of 
optical  scattering.  Turbidity  is  a  derived  quantity  which  can  also  be 
represented  as  an  empirical  function  of  suspended  solids  concentration. 
The  option  to  compute  turbidity  will  be  provided  in  the  GUI. 

3.5  Phytoplankton 

Algal  sources  and  sinks  in  the  conservation  equation  include  production, 
respiration,  predation,  and  settling.  These  are  expressed: 


o  o 

—B=  G-R-Wa  —  B-PR 
5f  5z 


(11) 
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in  which: 

B  =  algal  biomass,  expressed  as  carbon  (g  C  m-3) 

G  =  growth  (d-i) 

R  =  respiration  (d  ^) 

Wa  =  algal  settling  velocity  (m  d  ^) 

PR  =  predation  (g  C  m-3  d'O 
z  =  vertical  coordinate. 

3.5.1  Production 

Production  by  phytoplankton  is  determined  by  the  intensity  of  light,  by  the 
availability  of  nutrients,  and  by  the  ambient  temperature. 

3.5.2  Light 

The  influence  of  light  on  phytoplankton  production  is  represented  by  a 
chlorophyll-specific  production  equation  (Jassby  and  Platt  1976): 

P^  =  P^m-  ,  ^  (12) 


in  which: 

P^  =  photosynthetic  rate  (g  C  g-^  Chi  d'O 
P^m  =  maximum  photosynthetic  rate  (g  C  g-^  Chi  d'O 
I  =  irradiance  (E  m-^  d'O- 

Parameter  Ik  is  defined  as  the  irradiance  at  which  the  initial  slope  of  the 
production  vs.  irradiance  relationship  (Figure  5)  intersects  the  value  of 
PBjn: 


a 


(13) 


in  which: 


a  =  initial  slope  of  production  vs.  irradiance  relationship  (g  C  g-i 
Chi  (E  m-2)-i) 
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Figure  5.  The  production  vs.  irradiance  curve.  When  I  «  ik,  production  depends  on 
avaiiabie  iight.  When  i  »  Ik,  production  approaches  a  maximum  independent  of 

available  light. 


The  chlorophyll-specific  production  rate,  P^,  is  readily  converted  to  carbon 
specific  growth  rate,  for  use  in  Equation  ii,  through  division  by  the 
carbon-to-chlorophyll  ratio: 


pB 

CChl 


(14) 


in  which: 

CChl  =  carbon-to-chlorophyll  ratio  (g  C  g  ^  chlorophyll  a). 

3.5.3  Nutrients 

Carbon,  nitrogen,  and  phosphorus  are  the  primary  nutrients  required  for 
algal  growth.  Diatoms  require  silica,  as  well.  Inorganic  carbon  and  silica  are 
usually  available  in  excess  and  are  not  considered  in  the  model.  The  effects 
of  the  remaining  nutrients  on  growth  are  described  by  the  formulation 
commonly  referred  to  as  “Monod  kinetics”  (Figure  6;  Monod  1949): 
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fW 


D 

KHd  +  D 


(15) 


in  which: 

/(AO  =  nutrient  limitation  on  algal  production  (o  < /(AO  <  i) 

D  =  concentration  of  dissolved  nutrient  (g  m-3) 

KHd  =  half-saturation  constant  for  nutrient  uptake  (g  m-3). 

Figure  6.  Monod  formulation  for  nutrient-limited  growth.  N/KH  is  the  ratio  of  available 
nutrient  concentration  to  the  half-saturation  concentration  for  nutrient  uptake. 


N/KH 


3.5.4  Temperature 

Algal  production  increases  as  a  function  of  temperature  until  an  optimum 
temperature  or  temperature  range  is  reached.  Above  the  optimum, 
production  declines  until  a  temperature  lethal  to  the  organisms  is  attained, 
Numerous  functional  representations  of  temperature  effects  are  available. 
Inspection  of  growth  vs.  temperature  data  indicates  a  function  similar  to  a 
Gaussian  probability  curve  (Figure  7)  provides  a  good  fit  to  observations: 


f{T)  -  e^’^sMT-Toptf 

(16) 

f{T)  =  e-^92.(ropf-rf 

(17) 
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in  which: 

T  =  temperature  (oC) 

Topt  =  optimal  temperature  for  algal  growth  (oC) 

KTgi  =  effect  of  temperature  below  Topt  on  growth  (oC-2) 
KTg2  =  effect  of  temperature  above  Topt  on  growth  (oC  ^). 


Figure  7.  Dependence  of  algal  production  on  temperature. 


3.5.5  Combining  the  effects  of  iight,  nutrients,  and  temperature 

A  production  vs.  irradiance  relationship  (Equation  12)  is  constructed  for 
each  model  cell  at  each  time-step.  First,  the  maximum  photosynthetic  rate 
under  ambient  temperature  and  nutrient  concentrations  is  determined: 

P^m{N,T)  -  P^m ■  f(T) ■  f(N)  (18) 

in  which: 

P^m{N,T)  =  maximum  photosynthetic  rate  under  ambient  temperature  and 
nutrient  concentrations  (g  C  g  ^  Chi  d  ^) 

The  most  limiting  nutrient  is  employed  in  determining  the  nutrient 
limitation, /(AO-  Next,  parameter  Ik  is  derived  from  Equation  13.  Finally, 
the  production  vs.  irradiance  relationship  is  constructed  using  P^m{N,T) 
and  Ik.  The  production  vs.  irradiance  relationship  allows  the  calculation  of 
growth  rate,  G,  at  any  depth  and  level  of  irradiance  via  Equation  14. 
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3.5.6  Respiration 

Two  forms  of  respiration  are  considered  in  the  model:  photo-respiration 
and  basal  metabolism.  Photo-respiration  represents  the  energy  expended 
by  carbon  fixation  and  is  a  fixed  fraction  of  production.  In  the  event  of  no 
production  (e.g.,  at  night),  photo-respiration  is  zero.  Basal  metabolism  is  a 
continuous  energy  expenditure  to  maintain  basic  life  processes.  In  the 
model,  metabolism  is  considered  to  be  an  exponentially  increasing 
function  of  temperature  (Figure  8).  Total  respiration  is  represented: 

R  =  Presp  G  +  BM-  ^ ^  ( 19) 


in  which: 

Presp  =  photo-respiration  (o  <  Presp  <  i) 

BM  =  metabolic  rate  at  reference  temperature  Tr  (d  O 
KTb  =  effect  of  temperature  on  metabolism  (oC  O 
Tr  =  reference  temperature  for  metabolism  (oC). 


Figure  8.  Exponential  temperature  relationship  employed  for  algal  metabolism  and 

other  processes. 


T  -  Tref 


3.5.7  Predation 

Predation  is  modeled  by  assuming  predators  clear  a  specific  volume  of 
water  per  unit  biomass: 


PR  =  FBM 


(20) 
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in  which: 

F  =  filtration  rate  (m3  g  i  predator  C  d  ^) 

M  =  planktivore  biomass  (g  C  m-3). 

Detailed  specification  of  the  spatial  and  temporal  distribution  of  the 
predator  population  is  impossible.  One  approach  is  to  assume  predator 
biomass  is  proportional  to  algal  biomass,  M  =  y  5,  in  which  case  Equation 
20  can  be  rewritten: 


PR  =  yFB^  (21) 

Equation  2i  results  in  a  quadratic  predation  term;  predation  rate  is 
proportional  to  the  algal  biomass  squared.  This  formulation  is 
recommended  by  Cerco  and  Noel  (2004)  based  on  comparisons  of 
computed  and  observed  algal  biomass  and  production  in  Chesapeake  Bay. 
To  provide  maximum  model  flexibility,  however,  the  exponent  on  algal 
biomass  is  a  user-specified  variable.  The  predation  term  as  expressed  in  the 
model  code  follows: 


PR  =  BPRB^^^ 


(22) 


in  which: 

BPR  =  basic  predation  rate  (m3  g  i  C  d  ^) 

PRP  =  predation  power. 

The  value  PRP  =  2  results  in  the  recommended  formulation.  The  value 
PRP  =  1  results  in  a  first-order  loss  term  which  is  in  common  usage.  The 
units  of  BPR  vary,  depending  on  the  value  of  PRP.  The  units  specified 
above  are  for  PRP  =  2. 

3.5.8  Accounting  for  algal  nitrogen 

Model  nitrogen  state  variables  include  ammonium,  nitrate,  and  organic 
nitrogen.  The  amount  of  nitrogen  incorporated  in  algal  biomass  is 
quantified  through  a  stoichiometric  ratio.  Thus,  total  nitrogen  in  the 
model  is  expressed: 


TotN  =  NH^  +  ATOg  +  Anc  ■  B  +  OrgN 


(23) 
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in  which: 

TotN  =  total  nitrogen  (g  N  ni-3) 

NH4  =  ammonium  (g  N  m-3) 

ATOs  =  nitrate  (g  N  m-3) 

Anc  =  algal  nitrogen-to-carbon  ratio  (g  N  g-i  C) 
OrgN  =  organic  nitrogen  (g  N  m-3). 


Algae  take  up  ammonium  and  nitrate  during  production  and  release 
ammonium  and  organic  nitrogen  through  respiration.  The  fate  of  nitrogen 
released  by  respiration  is  determined  by  empirical  distribution 
coefficients.  The  fate  of  algal  nitrogen  recycled  by  predation  is  determined 
by  a  second  set  of  distribution  parameters. 


3.5.9  Algal  nitrogen  preference 


Algae  can  utilize  ammonium  and  nitrate  for  production.  Trace 
concentrations  of  ammonium  inhibit  nitrate  uptake  so  that,  in  the 
presence  of  multiple  nitrogenous  nutrients,  ammonium  is  utilized  first. 
The  preference  of  algae  for  ammonium  is  expressed  by  a  modification  of 
an  empirical  function  presented  by  Thomann  and  Fitzpatrick  (1982): 


PN 


=  NH, 


+NH,- 


_ ^ _ 

(KHNH^  +  NH^  )  ■  (KHNH^  + 

KHNH^ 

(NH^  +  ATOg )  ■  (KHNH^  + 


(24) 


in  which: 

PN  =  algal  preference  for  ammonium  uptake  (o  <  PN  <  1) 
KHNH4  =  half  saturation  concentration  for  algal  ammonium  uptake 
(g  N  m-3). 


The  preference  function  has  two  limiting  values  (Figure  9).  When  nitrate  is 
absent,  the  preference  for  ammonium  is  unity.  When  ammonium  is  absent, 
the  preference  is  zero.  In  the  presence  of  ammonium  and  nitrate,  the 
preference  depends  on  the  abundance  of  both  forms  relative  to  the  half¬ 
saturation  constant  for  nitrogen  uptake.  When  both  ammonium  and  nitrate 
are  abundant,  the  preference  for  ammonium  approaches  unity.  When 
ammonium  is  scarce  but  nitrate  is  abundant,  the  preference  decreases  in 
magnitude,  and  a  significant  fraction  of  algal  nitrogen  requirement  comes 
from  nitrate. 
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Figure  9.  Algal  ammonium  preference.  The  preference  depends  on  the  abudance  of 
ammonium  and  nitrate  relative  to  the  half-saturation  concentration  for  algal 


3.5.10  Accounting  for  algal  phosphorus 

As  with  nitrogen,  the  amount  of  phosphorus  incorporated  in  algal  biomass 
is  quantified  through  a  stoichiometric  ratio.  Thus,  total  phosphorus  in  the 
model  is  expressed: 


TotP  =  PO^  +  Ape  ■  B  +  OrgP 


(25) 


in  which: 

TotP  = 
PO4  = 
Ape  = 
OrgP  = 

total  phosphorus  (g  P  m-3) 
dissolved  phosphate  (g  P  m-3) 
algal  phosphorus-to-carbon  ratio  (g  P  g-^  C) 
organic  phosphorus  (g  P  m-3). 

Algae  take  up  phosphate  during  production  and  release  phosphate  and 
organic  phosphorus  through  respiration.  The  fate  of  phosphorus  released 
by  respiration  is  determined  by  empirical  distribution  coefficients.  The 
fate  of  algal  phosphorus  recycled  by  predation  is  determined  by  a  second 
set  of  distribution  parameters. 

3.5.11  Effect  of  algae  on  dissolved  oxygen 

Algae  produce  oxygen  during  photosynthesis  and  consume  oxygen 
through  respiration.  The  quantity  produced  depends  on  the  form  of 
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nitrogen  utilized  for  growth.  When  ammonium  is  the  nitrogen  source,  one 
mole  oxygen  is  produced  per  mole  carbon  dioxide  fixed.  When  nitrate  is 
the  nitrogen  source,  1.3  moles  oxygen  are  produced  per  mole  carbon 
dioxide  fixed. 


In  addition  to  algal  respiration,  the  model  allows  for  a  portion  of  algal 
biomass  lost  to  predation  to  be  represented  as  oxygen  consumption.  This 
consumption  represents  respiration  by  predators.  Oxygen  consumption  is 
limited  to  the  amount  available.  In  the  event  dissolved  oxygen  (DO) 
concentration  is  insufficient  to  supply  the  respiratory  demand,  a  portion  of 
respiration  is  instead  allocated  to  release  of  organic  carbon.  The  equation 
that  describes  the  effect  of  algae  on  dissolved  oxygen  in  the  model  follows: 


—DO  = 

8f 

1.3-0.3- PN)G-(1-FC)R\A0CRB 
-FDOP  (l-FC)-  AOCR  ■  PR 


(26) 


in  which: 

FC  =  fraction  of  oxygen  consumption  recycled  as  organic  carbon 
due  to  low  DO  concentration  (o  <  FC  <  1) 

AOCR  =  dissolved  oxygen-to-carbon  ratio  in  respiration  (2.67  g  O2  g'^ 
C) 

FDOP  =  fraction  of  predation  expressed  as  oxygen  consumption  (o  < 
FDOP  <  1). 

The  quantity  (1.3  -  0.3  •  PN)  is  the  photosynthesis  ratio  and  expresses  the 
molar  quantity  of  oxygen  produced  per  mole  carbon  fixed.  The  photo¬ 
synthesis  ratio  approaches  unity  as  the  algal  preference  for  ammonium 
approaches  unity.  The  fraction  of  oxygen  consumption  recycled  as  organic 
carbon  is  represented: 


FC 


KHr 

KHr  +  DO 


(27) 


in  which: 


KHr  =  dissolved  oxygen  concentration  at  which  oxygen  consumption 
is  halved  (g  DO  m-3). 
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3.6  Organic  carbon 

Organic  carbon  undergoes  numerous  transformations  in  the  water 
column.  The  model  carbon  cycle  (Figure  lo)  consists  of  the  following 
elements: 


•  phytoplankton  production  and  excretion 

•  predation  on  phytoplankton 

•  heterotrophic  respiration 

•  settling. 


Figure  10.  The  model  carbon  cycle.  Phytoplankton  and  organic  carbon  are  model  state 

variables. 


dissolved 

inorganiccarbon 


external  loads 


settling 


Algal  production  is  frequently  the  primary  carbon  source  although  carbon 
also  enters  the  system  through  external  loading.  Predation  on  algae  by 
zooplankton  and  other  organisms  releases  organic  carbon  to  the  water 
column.  Organic  carbon  is  respired  at  a  first-order  rate  to  inorganic 
carbon.  A  portion  of  the  organic  carbon  which  does  not  undergo 
respiration  settles  to  the  bottom  sediments. 

Organic  carbon  respiration  is  represented  as  a  first-order  process  in  which 
the  reaction  rate  is  proportional  to  concentration  of  the  reactant.  The 
reaction  is  restricted  at  low  DO  concentrations.  An  exponential  function 
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(Figure  8)  relates  respiration  to  temperature.  The  complete  representation 
of  organic  carbon  sources  and  sinks  in  the  model  system  follows: 


in  which: 


OrgC  =  organic  carbon  (g  m-3) 

FC  =  fraction  of  algal  and  predator  respiration  released  as  OrgC 
due  to  low  DO  concentration  (o  <  FC  <  i) 

R  =  algal  respiration  rate  (d  O 
B  =  algal  biomass  (g  m-3) 

FDOP  =  fraction  of  predation  expressed  as  oxygen  consumption  (o  < 
FDOP  <  1) 

PR  =  predation  on  algae  (g  C  m-3  d-^) 

KHr  =  dissolved  oxygen  concentration  at  which  respiration  of 
organic  carbon  is  halved  (g  DO  m-3) 

Korgc  =  respiration  rate  of  OrgC  (d-i) 

Wop  =  settling  rate  for  organic  matter  (m  d-i) 
z  =  vertical  coordinate. 


Predation  contributes  to  OrgC  through  two  terms.  The  first  represents  the 
fraction  of  predation  intended  for  oxygen  consumption,  but  for  which,  the 
dissolved  oxygen  supply  is  insufficient.  The  second  term  is  the  fraction 
routed  directly  to  OrgC  under  all  conditions. 

3.7  Nitrogen 

The  model  nitrogen  cycle  (Figure  ii)  includes  the  following  processes: 

•  algal  production  and  metabolism 

•  predation 

•  mineralization  of  organic  nitrogen 

•  settling 

•  nitrification. 
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Figure  11.  The  model  nitrogen  cycle.  Nitrate,  ammonium,  organic  nitrogen  and  phytoplankton 

are  model  state  variables. 


denitrification  sediment 

in  sediments  release 


settling 


External  loads  are  the  ultimate  source  of  nitrogen  to  the  system.  Available 
nitrogen  (ammonium  and  nitrate)  is  incorporated  by  algae  during  growth 
and  released  as  ammonium  and  organic  nitrogen  through  respiration  and 
predation.  A  portion  of  the  organic  nitrogen  mineralizes  to  ammonium.  An 
additional  portion  settles  to  the  sediments.  In  an  oxygenated  water 
column,  a  fraction  of  the  ammonium  is  subsequently  oxidized  to  nitrate 
through  the  nitrification  process.  Organic  nitrogen  that  settles  to  the 
sediments  is  mineralized  and  recycled  to  the  water  column,  primarily  as 
ammonium.  Nitrate  moves  in  both  directions  across  the  sediment-water 
interface,  depending  on  relative  concentrations  in  the  water  column  and 
sediment  interstices. 

3.7.1  Nitrification 

Nitrification  is  a  process  mediated  by  specialized  groups  of  autotrophic 
bacteria  that  obtain  energy  through  the  oxidation  of  ammonium  to  nitrite 
and  oxidation  of  nitrite  to  nitrate.  A  simplified  expression  for  complete 
nitrification  (Tchobanoglous  and  Schroeder  1987)  follows: 

yields 

NH+  +  2O2  ^  NO^  +  H^O  +  2H+  (29) 

The  simplified  stoichiometry  indicates  that  two  moles  of  oxygen  are 
required  to  nitrify  one  mole  of  ammonium  into  nitrate.  The  simplified 
equation  is  not  strictly  true,  however.  Cell  synthesis  by  nitrifying  bacteria 
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is  accomplished  by  the  fixation  of  carbon  dioxide  so  that  less  than  two 
moles  of  oxygen  are  consumed  per  mole  ammonium  utilized  (Wezernak 
and  Gannon  1968). 


The  kinetics  of  complete  nitrification  are  modeled  as  a  function  of 
available  ammonium,  dissolved  oxygen,  and  temperature: 


NT 


DO  NH^ 

KHont  +  DO  KHnnt  +  NH^ 


■f(T)NTm 


(30) 


in  which: 

NT  =  nitrification  rate  (g  N  m-3  d  ^) 

KHont  =  half-saturation  constant  of  dissolved  oxygen  required  for 
nitrification  (g  O2  m-3) 

KHnnt  =  half-saturation  constant  of  NH4  required  for  nitrification 
(g  N  m-3) 

NTm  =  maximum  nitrification  rate  at  optimal  temperature 
(g  N  m-3  d-i). 


The  kinetics  formulation  (Figure  12)  incorporates  the  products  of  two 
Monod-like  functions.  The  first  function  diminishes  nitrification  at  low 
dissolved  oxygen  concentration.  The  second  function  expresses  the 
influence  of  ammonium  concentration  on  nitrification.  When  ammonium 
concentration  is  low,  relative  to  KHnnt,  nitrification  is  proportional  to 
ammonium  concentration.  For  AT/q  <<  KHnnt,  the  reaction  is 
approximately  first-order.  (The  first-order  decay  constant!^  NTm/ KHnnt.) 
When  ammonium  concentration  is  large,  relative  to  KHnnt,  nitrification 
approaches  a  maximum  rate.  This  formulation  is  based  on  a  concept 
proposed  by  Tuffey  et  al.  (1974).  Nitrifying  bacteria  adhere  to  benthic  or 
suspended  sediments.  When  ammonium  is  scarce,  vacant  surfaces  suitable 
for  nitrifying  bacteria  exist.  As  ammonium  concentration  increases, 
bacterial  biomass  increases,  vacant  surfaces  are  occupied,  and  the  rate  of 
nitrification  increases.  The  bacterial  population  attains  maximum  density 
when  all  surfaces  suitable  for  bacteria  are  occupied.  At  this  point, 
nitrification  proceeds  at  a  maximum  rate  independent  of  additional 
increase  in  ammonium  concentration. 
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Figure  12.  Effect  of  dissolved  oxygen  and  ammonium  concentrations  on  nitrification 
rate.  When  ammonium  is  abundant,  NH4  »  KHnt,  nitrification  rate,  NT,  approaches 
a  maximum  value.  As  dissolved  oxygen  diminishes,  DO  «  KHnt,  nitrification  rate 

approaches  zero. 
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The  optimal  temperature  for  nitrification  may  be  less  than  peak 
temperatures  that  occur  in  surface  waters.  To  allow  for  a  decrease  in 
nitrification  at  superoptimal  temperature,  the  effect  of  temperature  on 
nitrification  is  modeled  in  the  Gaussian  form  of  Equations  16  and  17. 

3.7.2  Nitrogen  mass  balance  equations 

The  mass-balance  equations  for  nitrogen  state  variables  are  written  by 
summing  all  previouslydescribed  sources  and  sinks: 

3.7. 2.1  Ammonium 


—NH,= 

8f  " 

Anc  [(R- FNI  -  PN  G)- B  +  PR- FNIP 
■^Korgn  ■  OrgN  —  NT 


(31) 


in  which: 

Anc  =  algal  nitrogen-to-carbon  ratio  (g  N  g  i  C) 

R  =  algal  respiration  rate  (d-i) 

FNI  =  fraction  of  algal  metabolism  released  as  NH4  (o  <  FNI  <  1) 
PN  =  algal  ammonium  preference  (o  <  PN  <  1) 
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G 

B 

PR 

FNIP 

Korgn 

OrgN 


algal  growth  rate  (d-i) 

algal  biomass  (g  m-3) 

predation  on  algae  (g  C  m-3  d'^) 

fraction  of  predation  released  as  NH4  (o  <  FNIP  <  1) 

organic  nitrogen  mineralization  rate  (d'O 

organic  nitrogen  (g  m-3). 


3.7. 2.2  Nitrate 


^NO^=-Anc-(l-PN)-PB  +  NT 


(32) 


in  which: 


NO3  =  nitrate  concentration  (g  m-3). 
3. 7.2.3  Organic  Nitrogen 


8f 


OrgN 


Anc[RB(l-FNI)  +  PR(l-FNIP)]-KorgnOrgN 

o 

—Wop — OrgN 
5z 


(33) 


in  which: 

Wop  =  settling  rate  for  organic  matter  (m  d-^) 
z  =  vertical  coordinate. 

3.8  Phosphorus 

The  model  phosphorus  cycle  (Figure  13)  includes  the  following  processes: 

•  algal  uptake  and  excretion 

•  predation 

•  mineralization  of  organic  phosphorus 

•  settling. 
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Figure  13.  The  model  phosphorus  cycle.  Phosphate,  organic  phosphorus,  and  phytoplankton 

are  model  state  variables. 
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External  loads  are  the  ultimate  source  of  phosphorus  to  the  system. 
Dissolved  phosphate  is  incorporated  by  algae  during  growth  and  released, 
along  with  organic  phosphorus,  through  respiration  and  predation.  A 
portion  of  the  organic  phosphorus  is  mineralized  to  phosphate.  A  second 
portion  settles  to  the  sediments.  Within  the  sediments,  organic  phosphorus 
is  mineralized  and  recycled  to  the  water  column  as  dissolved  phosphate. 

3.8.1  Phosphate 

The  mass  balance  equation  for  phosphate,  including  all  processes,  is 

o 

— PO4  =  Korgp  ■  OrgP  -  Ape  - G  B 

5f  "  y  p  ^34^ 

+Apc  ■  [FBI  ■RB  +  FPIP  ■  PR] 


in  which: 


PO4 

Korgp 

OrgP 

Ape 


G 

B 

FPI 


R 

FPIP 


phosphate  (g  m-3) 

organic  phosphorus  mineralization  rate  (d  O 
organic  phosphorus  (g  m-3) 
algal  phosphorus-to-carbon  ratio  (g  P  g  ^  C) 
algal  growth  rate  (d  ^) 
algal  biomass  (g  m-3) 

fraction  of  algal  metabolism  released  as  phosphate 
(o  <PP/<  1) 

algal  respiration  rate  (d-^) 

fraction  of  predation  released  as  phosphate  (o  <  FPIP  <  1) 
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PR  =  predation  on  algae  (g  C  ni-3  d-i)- 

3.8.2  Organic  phosphorus 

The  mass  balance  equation  for  organic  phosphorus,  including  all 
processes, is 

o 

—OrgP  =  Apc[RB(l-FPI)  +  PR(l-FPIP)] 

(35) 

—Korgp  ■  OrgP  —  Wop  — OrgP 

8z 

in  which: 

Wop  =  settling  rate  for  organic  matter  (m  d  ^ 
z  =  vertical  coordinate. 

3.9  Dissolved  oxygen 

Sources  and  sinks  of  dissolved  oxygen  in  the  water  column  (Figure  14) 
include 

•  algal  photosynthesis 

•  atmospheric  reaeration 

•  algal  respiration 

•  heterotrophic  respiration 

•  nitrification. 
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Figure  14.  The  model  dissolved  oxygen  cycle.  Phytoplankton,  dissolved  oxygen,  ammonium, 
and  organic  carbon  are  model  state  variables. 
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3.9.1  Reaeration 

The  rate  of  reaeration  is  proportional  to  the  dissolved  oxygen  deficit  in 
model  segments  which  form  the  air-water  interface: 

8  Kr 

—DO  =  —(DOs-DO)  (36) 

8f  Az 


in  which: 

DO  =  dissolved  oxygen  concentration  (g  O2  m-3) 

Kr  =  reaeration  coefficient  (m  d  ^) 

DOs  =  dissolved  oxygen  saturation  concentration  (g  O2  m-3) 

Az  =  model  surface  layer  thickness  (m). 

In  free-flowing  streams,  the  reaeration  coefficient  depends  largely  on 
turbulence  generated  by  bottom  shear  stress  (O'Connor  and  Dobbins 
1958).  In  lakes  and  coastal  waters,  however,  wind  effects  dominate  the 
reaeration  process  (O'Connor  1983).  For  the  initial  version  of  ICM-Lite,  a 
relationship  for  wind-driven  gas  exchange  (Hartman  and  Hammond  1985) 
is  employed: 


Kr  =  Arear  ■  Rnu  ■  Wms^'^ 


(37) 


in  which: 
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Arear  =  empirical  constant  («  o.i) 

Rnu  =  ratio  of  kinematic  viscosity  of  pure  water  at  20  oC  to  kinematic 
viscosity  of  water  at  specified  temperature  and  salinity 
Wms  =  wind  speed  measured  at  lo  m  above  water  surface  (m  s  ^)- 

Hartman  and  Hammond  (1985)  indicate  Arear  takes  the  value  0.157.  In 
the  present  model,  Arear  is  treated  as  a  variable  to  allow  for  effects  of 
wind  sheltering,  for  differences  in  height  of  local  wind  observations,  and 
for  other  factors.  An  empirical  function  that  fits  tabulated  values  of  Rnu  is 

Rnu  =  0.54 +  0.0233-r-0.002-S  (38) 


in  which: 

S  =  salinity  (ppt) 

T  =  temperature  (oC). 

Saturation  dissolved  oxygen  concentration  diminishes  as  temperature  and 
salinity  increase.  An  empirical  formula  that  describes  these  effects  (Genet 
et  al.  1974)  is 


DOs  =  14.5532  -  0.382 17  •  T  +  0.0054258  • 

-CL  ■  (^1.665  X 10^^  -  5.866  x  10  ®  •  T  +  9.796  x  10  ®  •  T^) 


(39) 


in  which: 

CL  =  chloride  concentration  (=  salinity/1.80655). 

3.9.2  Mass  balance  equation  for  dissolved  oxygen 

The  summary  of  all  terms  which  affect  dissolved  oxygen  is 

o 

-^DO  =  AOCR  [{1.3 -0.3- PN)  G-(1-FC)  R\  B 

-FDOP  (l-FC)-  AOCR  ■  PR 

-AONT  ■  NT - — - AOCR  ■  Korgc  ■  OrgC 

KHr  +  DO 

+—(DOs-DO) 

Az 


(40) 


in  which: 
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AOCR  =  oxygen-to-carbon  mass  ratio  in  production  and  respiration 
(=  2.67  g02  g-i  C) 

PN  =  algal  ammonium  preference  (o  <  PN  <  1) 

G  =  algal  growth  rate  (d-i) 

FC  =  fraction  of  algal  and  predator  respiration  released  as  OrgC 
due  to  low  DO  concentration  (o  <  FC  <  1) 

R  =  algal  respiration  rate  (d-i) 

B  =  algal  biomass  (g  m-3) 

FDOP  =  fraction  of  predation  expressed  as  oxygen  consumption 
(o  <  FDOP  <  1) 

PR  =  predation  on  algae  (g  C  m-3  d'^) 

AONT  =  oxygen  consumed  per  mass  ammonium  nitrified 
(=  4-33  g  O2  g-i  N) 

NT  =  nitrification  rate  (g  N  m-3  d-^) 

KHr  =  dissolved  oxygen  concentration  at  which  respiration  of 
organic  carbon  is  halved  (g  DO  m-3) 

Korgc  =  respiration  rate  of  OrgC  (d-^) 

OrgC  =  organic  carbon  (g  m-3). 

3.10  User-defined  substance 

ICM-Lite  provides  a  user-defined  substance.  Basic  kinetics  including  first- 
order  decay  and  settling  are  considered: 

o  o 

—UDS  =  -Kuds  ■  UDS  -  Wuds—UDS  (41) 

5f  5z 


in  which: 

UDS  =  user-defined  substance  (mass  per  unit  volume) 
Kuds  =  decay  rate  (d-^) 

Wuds  =  settling  rate  (m  d-i)- 


The  concentration  of  UDS  will  usually  be  g  m-3,  consistent  with  most  other 
model  variables.  The  user  can  employ  alternate  units  while  taking  care  to 
employ  consistent  units  for  loads  and  boundary  conditions. 
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4  Sediment-Water  Interactions 

Exchange  of  material  between  the  water  column  and  benthic  sediments  is 
an  important  component  of  aquatic  ecosystems.  Over  lengthy  time  scales 
(e.g.,  years  to  decades),  the  sediments  are  an  ultimate  sink  of  nutrients 
and  organic  matter  which  settle  from  the  water  column.  Over  lesser  time 
scales  (e.g.,  seasons  to  years),  however,  sediment  release  of  previously 
deposited  nutrients  can  be  a  net  source  to  the  water  column.  Sediment 
oxygen  demand  may  comprise  a  substantial  fraction  of  total  system  oxygen 
consumption. 

The  complete  CE-QUAL-ICM  is  usually  executed  in  conjunction  with  a 
predictive  sediment  diagenesis  model  (DiToro  2001).  The  diagenesis  model 
is  demanding  in  terms  of  data  requirements  and  is  best  implemented  by  the 
experienced  user.  As  an  alternative,  ICM-Lite  implements  user-specified 
sediment-water  nutrient  and  oxygen  fluxes.  Basic  relationships  are  provided 
that  express  the  influence  of  conditions  in  the  water  column  on  the  specified 
fluxes. 

User-specified  fluxes  are  available  for 

•  ammonium 

•  nitrate 

•  phosphate 

•  dissolved  oxygen. 

Transfer  of  particulate  matter  from  water  to  sediments  is  treated  through 
specification  of  a  settling  velocity.  The  model  employs  the  convention  that 
positive  fluxes  are  from  sediment  to  water  and  negative  fluxes  are  from 
water  to  sediment.  Eluxes  of  ammonium  and  phosphate  are  most  often 
from  sediments  to  water  and  are  positive  quantities.  Nitrate  commonly 
passes  in  both  directions  across  the  sediment-water  interface  and  may  be 
positive  or  negative.  Since  oxygen  moves  from  water  to  sediments, 
sediment  oxygen  consumption  is  represented  as  a  negative  quantity. 


ERDC/EL  TR-15-8 


38 


4.1  Ammonium  and  phosphate 

Sediment  releases  of  ammonium  and  phosphate  are  specified  at  a  reference 
temperature.  The  effect  of  temperature  on  release  rate  (Figure  15)  is 
determined  by  an  exponential  relationship: 

BEN  =  BENb  •  (42) 


in  which: 

BEN  =  sediment  release  at  temperature  T  (g  m-2  d  O 
BENb  =  sediment  release  (g  m-2  d  ^)  specified  at  reference  temperature 
TRben 

KS  =  effect  of  temperature  on  sediment  release  (C°  ^ 

T  =  temperature  (0°). 


Figure  15.  Effect  of  temperature  on  sediment-water  fiuxes,  for  three  vaiues  of  KS.  When  T  = 
TRben,  the  fiux  equais  the  specified  base  vaiue. 


4.2  Nitrate 

Movement  of  nitrate  between  water  and  sediments  is  strongly  influenced 
by  concentration  of  nitrate  in  the  water  column.  When  nitrate  is  abundant 
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in  the  water  column,  nitrate  usually  diffuses  from  overlying  water  into  the 
sediments  where  it  is  denitrified  to  gaseous  form.  When  nitrate  is  absent 
from  the  water  column,  small  quantities  of  nitrate  may  diffuse  from 
sediment  interstitial  water  into  the  overlying  water.  The  model  allows  for 
user-specified  nitrate  flux  and  provides  a  function  (Figure  i6)  that  relates 
flux  to  concentration: 


BENNO^  =  BENNO^b 
+MTC  ■  (SEDNO^  -  iVOg  )  ■ 


(43) 


in  which: 


BENNO3  = 
BENN03b= 
MTC  = 
SEDNO3  = 
NO^  = 


sediment-water  nitrate  flux  (g  N  m-2  d  ^ 
specified  sediment-water  nitrate  flux  (g  N  m-2  d  O 
sediment-water  mass  transfer  coefficient  (m  d  ^) 
nitrate  concentration  in  interstitial  water  (g  m-3) 
nitrate  concentration  in  water  overlying  sediments  (g  m-3) 


Figure  16.  Effect  of  nitrate  concentration  in  water  column  on  sediment-water  nitrate  fiux. 
Caicuiated  for  MTC  =  0.1  m  d  1,  SEDN03  =  0.02  g  m  3.  When  NO3  =  SEDNO3,  the  fiux  is  zero. 
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In  model  employment,  the  user  specifies  the  nitrate  flux,  if  known,  or 
relies  on  the  model  to  compute  flux  as  a  function  of  nitrate  and 
temperature  in  the  water  column. 

4.3  Sediment  oxygen  consumption 

Oxygen  consumption  in  the  sediments  depends  upon  water-column 
temperature  and  oxygen  availability.  As  temperature  increases,  respiration 
in  the  sediment  increases.  Sediment  oxygen  consumption  is  reduced  as 
oxygen  concentration  in  the  overlying  water  decreases  (Figure  17).  The 
model  accounts  for  these  influences  through  the  relationship 

BENDO  = - — - BENDOb  ■  e®  (44) 

KHso  +  DO 

in  which: 

BENDO  =  sediment  oxygen  consumption  (g  m-2  d  ^) 

BENDOb  =  sediment  oxygen  consumption  under  conditions  of  unlimited 
oxygen  availability  (g  m-2  d-i),  specified  at  reference 
temperature  TRben 

KHso  =  dissolved  oxygen  concentration  at  which  sediment  oxygen 
consumption  is  halved  (g  m-3). 

4.4  Parameter  evaluation 

Base  fluxes  and  influences  of  temperature  and  other  factors  are  best 
determined  from  observations  collected  in  the  prototype  system.  Table  1 
lists  observations  from  several  systems  that  may  be  employed  as  starting 
values  when  no  observations  are  available.  Suggested  starting  values  for 
parameters  in  the  functions  that  relate  sediment-water  fluxes  to 
conditions  in  the  water  column  are  listed  in  Table  2. 
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Figure  17.  Effect  of  DO  concentration  on  sediment  oxygen  consumption.  Calcuiated  for  KHso 
=  1  g  m  3.  When  water  coiumn  dissoived  oxygen  equais  KHso,  the  sediment  oxygen 

consumption  is  haived. 


Tabie  1.  Observed  sediment-water  fiuxes. 


Ammonium, 

g  m-2  d'^ 

Nitrate, 

g  m-2  d'^ 

Phosphate, 

g  m-2  d'^ 

SOD, 

g  m-2  d'^ 

System 

0.01  to  0.28 

-0.04  to  0.1 

-0.003  to  0.03 

-1.5  to 
-3.5 

Chesapeake  Bay 
(Boynton  and 

Kemp  1985) 

-0.001  to  0.09 

-0.02  to  0.015 

-0.007  to  0.031 

-0.1  to 
-2.6 

Narragansett  Bay 
(Hale  1975) 

0  to  0.15 

0  to  0.002 

-0.006  to  0.034 

-0.6  to 
-2.4 

Neuse  and  South 
Rivers,  NO  (Fisher 
et  al.  1982) 

-0.04  to  0.36 

-0.1  to  0.08 

-0.019  to  0.124 

-0.1  to 
-2.7 

Potomac  Estuary 
(Callender  and 
Hammond  1982) 

-0.035  to  0.53 

-0.23  to  0.03 

0.001  to  0.22 

-0.5  to 
-4.1 

Patuxent  Estuary 
(Boynton  et  al. 

1980) 

Table  2.  Parameters  in  sediment-water  flux 


relationships. 


Parameter 

Suggested  Range 

KS 

0.04  to  0.07  °C-i 

MTC 

0.05  to  0.15  m  d'l 

SEDNO3 

Oto  0.05  g  m-3 

KHso 

1  to  2  g  m-3 
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